Stochastic switching of delayed feedback suppresses oscillations in genetic regulatory systems

Delays and stochasticity have both served as crucially valuable ingredients in mathematical descriptions of control, physical and biological systems. In this work, we investigate how explicitly dynamical stochasticity in delays modulates the effect of delayed feedback. To do so, we consider a hybrid model where stochastic delays evolve by a continuous-time Markov chain, and between switching events, the system of interest evolves via a deterministic delay equation. Our main contribution is the calculation of an effective delay equation in the fast switching limit. This effective equation maintains the influence of all subsystem delays and cannot be replaced with a single effective delay. To illustrate the relevance of this calculation, we investigate a simple model of stochastically switching delayed feedback motivated by gene regulation. We show that sufficiently fast switching between two oscillatory subsystems can yield stable dynamics.


Introduction
Mathematical models with delays have served profoundly useful in capturing the behaviour of complex systems in biology [1][2][3][4][5][6][7], networks [8] and control [9]. One notable example is delayed negative feedback control of genetic networks, especially transcriptional feedback [10][11][12]. These systems, such as in nF-κB [13] or p53 [14] signalling, share the canonical set-up of some signal that auto-inhibits its own production with delayed feedback. Such delays typically arise from several molecular events that must occur in sequence [15]. The behaviour of these models for a fixed τ is rich but well understood: the amount of delay in feedback crucially determines stability or instability (often to oscillations) in the system [1,16].
At the scale of molecular machinery associated with genetic regulation, dynamics are also known to be richly stochastic [17,18], with inherent noise in the counts of individual molecules, but also disparate timescales of promotion or inhibitory factors binding and unbinding [19]. Putting these two pieces together, there is a natural interest in understanding the emergent dynamics in systems with delayed feedback and stochasticity in biological [20][21][22][23][24][25][26][27][28] and physical systems [29,30].
There have been many insightful investigations into systems with distributed delay. In these models, a delay is continuously drawn from some probability distribution thought to arise from stochasticity or uncertainty in the delay [31,32]. Distributed delays in negative feedback can produce interesting behaviours, including bimodality [33] and surprising stability changes [11]. However, it remains unclear to what extent it is possible to relate distributed delays to more explicit dynamical descriptions of the underlying process that governs them. Others have considered discrete-time models or so-called semi-discretized where the delays themselves switch at discrete times [34][35][36][37][38]. Such models may be appropriate for control systems but unrealistic for biological systems where the switching is driven by stochastically timed chemical events [39]. Instead, we return to very early models where the delays themselves follow a continuous-time Markov process [40][41][42]. While there are rigorous works investigating the long-time stability of such models [43][44][45][46], these arguments are primarily based on Lyapunov functions that are challenging to find for any specific model.
In this work, we investigate the dynamics of a model that stochastically switches between delays at exponential rates. That is, we consider a stochastic hybrid delay system. The delays evolve via a continuous-time Markov chain, and in between these transitions, the system evolves by a deterministic delay equation. Our primary contribution is the derivation of an effective quasi-steady-state delay system in the limit of fast stochastic switching. Perhaps surprisingly, we find that nonlinear systems do not converge to one with a single effective delay but retain the effects of all delays in the original subsystems. Using this computation, we explore the behaviour of a classical model of delayed negative feedback [11,[47][48][49]. With stochastic switching, we show that switching between two oscillatory subsystems can stabilize the system. Altogether, our results add clarity and intrigue to the picture of how stochastic switching and delays intertwine in biological systems, especially those containing negative feedback as seen in genetic regulation.

Simple model of delayed negative feedback
We begin with a simple model of delayed negative feedback with a single fixed delay τ and review the role delay plays in destabilizing the fixed point of a dynamical system. We are far from the first to consider such a model or its variants [11,12,48,49], but we include a brief investigation here for the sake of self-containment of our work. Thereafter we demonstrate that allowing for stochastic switching between distinct delays stabilizes the fixed point, even when all delays involved are past the Hopf bifurcation point of the non-switching system.
Let y(t) be a scalar field evolving according to dy dt ¼ I À gy À wf ðyðt À tÞÞ: ð2:1Þ Here, y(t) could represent the concentration of a protein that is constitutively produced at a rate I and inhibits its own production. The first order rate constant γ describes the natural degradation rate of the substance y. We take the weight w to describe the strength of the auto-inhibition based on some form of Michaelis-Menten kinetics. Therefore, f can be any function that is monotonically increasing and finite at infinity. For concreteness, we take f ðyÞ ¼ y n K n þ y n : Equation (2.1) describes a simple model for delayed negative feedback-the substance being produced inhibits its own production following a non-zero but finite time delay τ.

Hopf bifurcation
Setting the time derivative equal to zero in equation (2.1) yields an equilibrium solution y* satisfying gy Ã þ wf ðy Ã Þ ¼ I: ð2:2Þ To understand how τ destabilizes the equilibrium, we perform a linear stability analysis. Linearizing equation (2.1) around y* yields where u(t) ≡ y(t) − y*. This has the solution u(t) = K e λt with λ determined from the eigenvalue equation l ¼ Àg À wf 0 ðy Ã Þ e Àlt : ð2:3Þ In accordance with standard analysis of delay differential equations, we determine the necessary conditions for the emergence of a time-periodic solution via a Hopf bifurcation by setting λ = iω in equation (2.3). Equating real and imaginary parts gives the following conditions for a Hopf bifurcation: v ¼ wf 0 ðy Ã Þ sin ðvtÞ Àg ¼ wf 0 ðy Ã Þ cos ðvtÞ It is clear that these conditions cannot be satisfied in the absence of delay (τ = 0). Indeed, setting τ = 0 renders equation (2.1) a one-dimensional flow, which cannot have oscillations. On the other hand, for τ > 0, there exist (ω c , τ c ) satisfying conditions (2.4). When τ is increased past τ c , a pair of complex conjugate eigenvalues crosses the imaginary axis. Although this is not sufficient to guarantee the emergence of a stable periodic solution via a supercritical Hopf bifurcation for τ > τ c , the existence of stable oscillations beyond the Hopf bifurcation point can be verified numerically, with examples seen in figure 1. In the simulations shown, the parameters chosen are I = 10, K = 9.5, γ = 1, n = 4 and w = 9.5. These will be used elsewhere throughout the manuscript unless noted otherwise. For these parameter values, τ c ≈ 0.9. The two delays in figure 1 correspond to values of τ above and below this bifurcation, τ = 0.6 and τ = 1.2, respectively.
Hence, we have established the crucial role delay plays in the destabilization of equilibria corresponding to delayed negative feedback systems. We next demonstrate the curious result that allowing τ to randomly switch between two values-both past τ c -results in the stabilization of the equilibrium.

Stochastic switching of delays
In biological and biophysical models, delays often manifest as a coarse-grained description of several processes which, cumulatively, require time. Hence, many details underlying biophysical processes are overlooked in fixed-delay systems. To capture finer dynamics, biophysical models incorporate delays that vary according to a probability distribution or evolve via some prescribed dynamics. We emphasize the contrast of this model set-up to those such as in [26,28], where stochastic switching of the system occurs (due to inhibitor binding and unbinding) but the delay, if included, remains fixed.
Here, we implement stochasticity in delay equations by taking the delay to evolve in time according to a continuous-time discrete Markov process. Explicitly, consider a general autonomous delay differential equation of the form dy dt ¼ GðyðtÞ, yðt À tðtÞÞÞ: ð2:5Þ Here y(t) is the same scalar field and τ(t) varies in time stochastically between N [ N states. Let τ i denote the delay corresponding to the ith state so that τ ∈ {τ 1 , τ 2 , …, τ N } at any given time t, and let Q(τ i , t) denote the probability that τ = τ i at time t. The dynamics of τ are then completely characterized by the master equation where W i!j denotes the propensity of the transition τ i → τ j . Equations (2.5) and (2.6) together form a so-called stochastic hybrid system-a system wherein the state of the system evolves stochastically, but within each state evolves deterministically. Such systems are also called a piecewise deterministic Markov process (PDMP) [19,39,50].

Stochastic delayed negative feedback
We now implement stochastic switching in equation (2.1) by taking τ to evolve according to a two-state Markov process, The last equality follows from the fact that Q(τ 1 , t) + Q(τ 2 , t) = 1 for this model. When τ is fixed, equation (2.1) has a welldefined Hopf bifurcation point at τ = τ c . Indeed, for τ 1 , τ 2 > τ c substituted into equation (2.1), the system admits a limit cycle with amplitude and frequency determined by the corresponding delay. Numerically simulating equation (2.1) with switching given in equation (2.7) indeed shows dynamics wherein the solution jumps between limit cycles (see figure 2b). However, in many situations, transitions between states of a PDMP are fast relative to the other dynamics of the system. Simulations of the stochastic hybrid system with τ 1 , τ 2 > τ c and a ¼ a=1 and b ¼ b=1 for 0 , 1 ( 1 show that the system contracts to the equilibrium (see figure 3). This result is not intuitive.

General formulation
To understand what causes the stabilization of the fixed point with fast stochastic delay switching, we consider a microscopic model for delayed negative feedback. We begin with a delayed master equation [51] and then invoke a van Kampen expansion to derive a Smoluchowski equation for the probability density of the stochastic variable undergoing delayed negative feedback. We begin by writing the master equation analogue of equation (2.1) with τ undergoing a continuous-time discrete Markov process with N states: Let dðtÞ [ N be a random variable whose dynamics are governed by the reactions delineated by equation (2.1). Thus, d(t) could represent the number of translated protein molecules. Let P(m, t, τ) be the probability that d = m at time t and that the current delay value is τ. The dynamics of P(m, t, τ) can be written as The first two terms correspond to constitutive protein production and natural degradation, respectively. The transitions between distinct delay values are captured in the last two terms of equation ( royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 20: 20230059 one needs to have information about three-and four-point correlations. The result is an infinite hierarchy of equations. Following a common approach [52,53], we assume the probability distributions for the number of protein molecules at distinct times are independent so that P(m, t, τ j ; M, t − τ i , τ j ) = P(m, t, τ j )P(M, t − τ i , τ j ). We refer the reader to the classical work of Frank [54,55] or the more recent [56] for further discussion on this approximation.
Unfortunately, little analysis can be performed on equation (3.1). To make progress, we follow arguments similar to those in [57] and invoke the van Kampen expansion by setting x ; m=N and Taylor expanding terms in equation (3.1) to OðN À1 Þ. Here, N represents the maximum number of translated protein molecules possible, as determined by the physics of the system. We obtain the following Chapman-Kolmogorov equations, as known for stochastic delay systems [56,58]: where p i (x, t) ≡ p(x, t, τ i ). The terms inside the summation and integral evaluate to the expected value of the Hill function, giving Here, χ i is the value of the delayed variable ξ yielding the expected value of the Hill function. It is uniquely identified because the Hill function is injective. We point out that equation (3.3) is a function of three independent variables: (x, χ i , t). Abusing notation slightly, we now set p i (x, t) ≡ p(x, χ i , t).
To determine why fast switching between delays stabilizes the fixed point of equation (2.1), we next must consider the limit where transitions between discrete delays are fast.

Fast switching limit
To consider the fast switching limit and derive the effective equation governing the dynamics of delayed negative feedback, we re-scale all transition propensities as W i!j ! 1 À1 W i!j , with 0 , 1 ( 1. The effective equation manifests as a perturbation from the stationary measure of the Markov transition matrix governing the switching between delays. It will provide an approximation to the mean dynamics of equation (3.3). First, we rewrite equation (3.3) in matrix-vector format:  Figure 2. Delayed negative feedback model with stochastic switching. (a) Schematic for the biological motivation of the model. A protein with concentration y(t) inhibits its own production with delay that now stochastically switches between τ 1 , τ 2 at exponentially distributed times. (b) Illustrative simulation of the system, here with relatively slow switching 1 ¼ 10. In the red regions, the delay state is τ 1 = 3 and in the blue, τ 2 = 1.  royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 20: 20230059 The linear operators L i act upon any sufficiently smooth function F by The co-kernel of A is spanned by the N-dimensional vector hcj ; ð1, 1, . . . , 1Þ, and the kernel of A is spanned by the N-dimensional vector |ϕ〉 ≡ (ϕ 1 , ϕ 2 , …, ϕ N ) T , whose entries, in general, will be determined by specific relations between the transition propensities and satisfy 〈ψ|ϕ〉 = 1. Indeed, |ϕ〉 is the invariant measure of the continuous-time Markov chain governing the jumps between delays. Let q = 〈ψ|p〉 and |w〉 = |p〉 − q|ϕ〉 so that q is the component of |p〉 in the co-kernel of A and |w〉 is in the orthogonal complement. Applying 〈ψ| to both sides of equation (3.4) gives @q @t ¼ hcjLjw þ qfi: ð3:5Þ Substituting |p〉 = |w + qϕ〉 into equation (3.4) and using the fact that |ϕ〉 is in the kernel of A, we obtain where I N is the N × N identity matrix. Introduce the expansion and substitute into equation (3.6). Collecting Oð1 À1 Þ terms gives A|w 0 〉 = |0〉 ↠ |w 0 〉 = |0〉. Because we are only interested in the mean dynamics and not the fluctuations about the mean, we need not attempt to calculate the higher-order terms. We hope to explore this in future work. Substituting the leading-order term for |w〉 into equation (3.5) gives the Smoluchowski equation for the probability density qðx, x, tÞ for a protein undergoing delayed negative feedback at time t, In words, equation (3.7) says that the leading-order effective dynamics evolve via the weighted average of the different subsystems, where the weight is determined by the stationary distribution of the underlying Markov chain controlling switching. An important caveat appears: that such a stationary distribution exists. For all examples within the realm of models of genetic feedback, we anticipate this to be the case, but the breakdown of this assumption may be of future mathematical interest.
Another area of mathematical interest is that equation (3.7) describes the dynamics of qðx, x, t ) with a multivariate Smoluchowski equation, therefore assuming that the current state of the random variable is independent of its value at all previous times. It is a Markovian description of a non-Markovian process-the fact that this is a delayed negative feedback system necessarily renders it non-Markovian. We find that equation (3.7) is a good approximation to the non-Markovian process in the fast switching limit (see figure 3d). We take this as evidence that the independence assumption we made works well in the fast switching limit. We hope to quantify this approximation and precisely determine where it breaks down in future work.
One notable feature of the structure of the effective dynamics (3.7) is that the system cannot be described by a single effective delay. Rather, the leading-order dynamics evolve by the simultaneous influence of all delay subsystems in the fast switching limit.

Delayed negative feedback again
We now invoke equation (3.7) and apply it to the specific example discussed in §2.3. In this case, the resulting Chapman-Kolmogorov equations are À ap 1 ðx, tÞ þ bp 2 ðx, tÞ and @p 2 @t þ ap 1 ðx, tÞ À bp 2 ðx, tÞ: Here, 〈ψ| = (1, 1) and so that the resulting Smoluchowski equation is qðx, x, tÞ : ð4:1Þ Thus, the mean dynamics of the protein undergoing delayed negative feedback is given by the kinetic equation wy eff ðt À t 1 Þ n K n þ y eff ðt À t 1 Þ n À a b þ a wy eff ðt À t 2 Þ n K n þ y eff ðt À t 2 Þ n : ð4:2Þ The stability of this averaged two-delay system can surprisingly display behaviour distinct from either of the two subsystems, as we show in the next section.
where uðtÞ ; y eff ðtÞ À y Ã . As in §2.1, we invoke the ansatz u(t) = A e λt and determine λ from the auxiliary equation Setting λ = iω yields the following conditions for a Hopf bifurcation in equation ð4:3Þ In figure 4, we show the locus of Hopf bifurcation points in parameter space for equation (4.2) and compare it with the Hopf bifurcation points for the single delay equation given in equation (2.1). We can see that there are regions of parameter space wherein τ 1 and τ 2 are larger than the single delay critical Hopf value but the system continues to reach the fixed point. Hence, the fast switching between delays in the stochastic system causes the effective behaviour of the system to behave as if the feedback followed two distinct delay values simultaneously. The presence of multiple delays increased the range of delay values for which the fixed point was stable. When switching is not fast, then the increased stabilization of the fixed point is not observed. Although it is challenging to analytically investigate this scenario, numerical investigation via stochastic simulation is straightforward and can be seen in figure 5. Explicit stochastic simulations are performed by sampling the continuous-time Markov chain and solving the delay differential equation between these events. The system is simulated for t ∈ [0, 100], and over the window t ∈ [90, 100], the minimum and maximum values are taken, as presumed magnitudes of any oscillations after the transient portions have decayed. As waiting times are increased (1 large), the system spends enough time in each delay state so that the effective dynamics follow a single delay equation for the duration of time in that state. Periodic solutions corresponding to the delay of the state emerge. As the waiting time is decreased (1 small), the effects of the second delay term emerge and the stabilization of the fixed point is observed (see figures 3 and 5).

Conclusion
We summarize the main contributions of the manuscript as follows. Most generally, we have derived an effective delay equation in the limit of fast switching between subsystems with different delays that evolve via a continuous-time Markov chain. A priori, it is not clear whether the behaviour when rapidly switching between systems can be replaced by one effective delay. Here, we answer that possibility in the negative for nonlinear systems, similar to the semi-discretized case [37]. We used this result to investigate a classical model of delayed negative feedback with a new twist of stochastic switching between two delays. In our stochastic model, we showed that sufficiently fast stochastic switching between two delays stabilizes the system where each delay alone produces oscillations.
Our results sit within broader biological and mathematical contexts. First, we note the relation to the literature on distributed delay systems, especially in models of genetic feedback. The effective dynamics derived here arising in the fast switching limit (4.2) are exactly the form of distributed delay descriptions of genetic feedback considered elsewhere [23,31,59]. We have, therefore, provided further mechanistic motivation for the inclusion of these distributed delay systems. We show that a Hopf bifurcation in the total switching rate occurs, indicating that fast switching and slow have fundamentally different behaviour. This nuance in timescales of stochasticity does not exist in descriptions with distributed delays.
Stochastically switching delays add a new vignette to the broader theme of stochasticity in genetic feedback. Importantly, we consider stochasticity only in the delay to emphasize its impact on the behaviour of the system. This is in contrast with other studies where stochasticity is included in genetic feedback in other ways and new behaviours appear. For instance, molecule counts in the genetic machinery are low enough to justify exploring demographic fluctuations [21]. Demographic noise can destabilize oscillations [27], whereas distributed delays can sharpen them  royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 20: 20230059 [60]. It is, therefore, of future interest to investigate how stochastic switching of delays and demographic fluctuations intertwine. A natural starting point would be similar investigations for non-delay systems, [61], but we anticipate further challenges due to the hierarchy of multi-point correlations for stochastic delay systems as discussed in [54][55][56][57]. Other studies also include stochastic switching with delays arising from binding and unbinding of promoters [26,28] producing bursting behaviour. These models have rich analytical tractability but do not include different delays. It remains to be explored whether these calculations can be extended to different delays, as considered here. To further the biological relevance, it is also pressing to develop a more mechanistic explanation of how stochastic switching of delays may arise. One intriguing direction is the emergence of stochastic switching from dual delay feedback pathways, such as in NF-κB signalling [49,62]. Lastly, on the purely mathematical side, our model and analysis add to the mosaic of stochastic systems that behave fundamentally differently than their deterministic counterparts or subsystems. Specifically, our results first add to a long history along the theme of how noise can stabilize systems [11,31,[63][64][65]. Secondly, they provide another example of how stochastic switching can result in unexpected behaviour of stable solutions to the corresponding non-switching ODEs and PDEs [66][67][68][69]. Although it may seem restrictive that we consider a Markov chain model by which the delays evolve, an arbitrary delay distribution may be constructed via the theory of phase-type distributions [70]. We were unable to compute analytical results when switching was not fast, or even a next-order correction to the leadingorder behaviour. It is perhaps feasible to use a momentbased approach that others have used in stochastically switching delayed [36,[71][72][73] or other stochastic hybrid systems [74,75]. Alternatively, it may be feasible and interesting to investigate the opposite limit of the one considered, where the dynamics of the subsystems are fast relative to the switching, as seen in many other biological systems [76,77], and is perhaps readily handled by classical homogenization techniques [78]. Further, the dynamical systems behaviour of the derived effective (4.2) system could also be probed more thoroughly, perhaps using the methods of Du et al. [79] that compute normal forms and investigate higher co-dimension bifurcations.
Data accessibility. All Matlab code for stochastic simulations, the Hopf bifurcation analysis, and generation of figures in the manuscript can be found from the GitHub repository: https://github.com/ chris-miles/switching-delays. only τ 2 only τ 1 log 10 (ε) log 10 (ε) max ω _S(ω)_ y max/min Figure 5. Bifurcation structure of the stochastically switching delay system as a function of the overall switching rate 1. Each circle is a stochastic realization, with 500 total per parameter set. (a) The maximum and minimum values over a time window in stochastic simulation show that as switching gets slower, oscillations emerge. Solid lines represent the oscillation peaks for the single delay subsystems with corresponding delays. (b) The peak of the power spectrum max v jSðvÞj also demonstrates a subcritical Hopf bifurcation parametrized by the switching rate 1.